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Abstract 

The numerical simulation of compactons, solitary waves with compact support, is 
. characterized by the presence of spurious phenomena, as numerically-induced radia- 

\ tion, which is illustrated here using four numerical methods applied to the Rosenau- 

Hyman K{p,p) equation. Both forward and backward radiations are emitted from 
the compacton presenting a self-similar shape which has been illustrated graphically 
by the proper scaling. A grid refinement study shows that the amplitude of the ra- 
^> ' diations decreases as the grid size does, confirming its numerical origin. The front 

■ velocity and the amplitude of both radiations have been studied as a function of 

T^lj- , both the compacton and the numerical parameters. The amplitude of the radiations 

I decreases exponentially in time, being characterized by a nearly constant scaling ex- 

OO ' ponent. An ansatz for both the backward and forward radiations corresponding to 

, a self-similar function characterized by the scaling exponent is suggested by the 

^ I present numerical results. 
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1 Introduction 



Compactons are travelling wave solutions with compact support resulting from 
the balance of both nonlinearity and nonlinear dispersion. Compacton solu- 
tions have been first found in a generalized Korteweg-de Vries equation with 
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nonlinear dispersion, the so-called (focusing) K{p,p) compacton equation of 
Rosenau and Hyman [1], given by 

du du du^ d^vP ^ 
dt dx dx dx^ 



where u{x^ t) is the wave amplitude, x is the spatial coordinate, t is time, 
and Co is a constant velocity used here in order to stop the compacton when 
required. Compactons are classical solutions of this equation only for 1 < p < 
3, otherwise they are "non-classical" or weak solutions. In fact, compacton 
solutions of Eq. (1), for p ^ {—1, 0, 1}, can be written as [2] 
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where c is the compacton velocity and Xq the position of its maximum at t = 0. 
Note that the compacton has k continuous derivatives at its both edges when 
p^{2 + k)/k. 

Compactons have multiple applications in Physics. Rosenau-Hyman (RH) 
equation (1) was discovered as a simplified model to study the role of nonlin- 
ear dispersion on pattern formation in liquid drops [1], being also proposed in 
the analysis of patterns on liquid surfaces [3]. Equations with compacton so- 
lutions have also found applications such as the lubrication approximation for 
thin viscous films [4], semiclassical models for Bose- Einstein condensates [5], 
long nonlinear surface waves in a rotating ocean when the high-frequency 
dispersion is null [6], the pulse propagation in ventricle-aorta system [7], dis- 
persive models for magma dynamics [8], or, even, particle wavefunctions in 
nonlinear quantum mechanics [9]. RH equation is also the continuous limit 
of the discrete equations of a nonlinear lattice [1]. In nonlinear lattices the 
propagation of compacton-like kinks has been observed using mechanical [10], 
electrical [11,12], and magnetic [13] analogs. Recently, RH equation has been 
generalized using a cosine nonlinearity in order to model the dispersive cou- 
pling in chains of oscillators resulting in the so-called phase compactons and 
kovatons [14], having apphcation in superconducting Josephson junction trans- 
mission lines [15]. Finally, let us remark that the general K{p, q) equation, with 
p ^ q, may also show elliptic function compactons [1,16,17], and that recent 
interest is focusing on multidimensional compactons [18,19]. 

The numerical simulation of the propagation of nonlinear waves presents sev- 
eral numerically induced phenomena, such as spurious radiation, artificial dis- 
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sipation, and errors in group velocity. The numerical analysis of compactons 
are not free of these spurious phenomena. In fact, the numerical solution of 
compacton equations is a very challenging problem presenting several numer- 
ical difficulties which has not been currently explained [16,20,21]. 

The numerical simulation of compactons by means of pseudospectral methods 
in space require the addition of artificial dissipation (hyperviscosity) using 
high-pass filters [1,19,22,23] in order to obtain stable results without appre- 
ciable spurious radiation. In fact, using those methods, Ref. [1] shows that 
compactons collide elastically, without visible radiation. However, after the 
collision, compactons show a phase shift and a small-amplitude, zero-mass, 
compact ripple is generated, which slowly decomposes into tiny compacton- 
anticompacton pairs. Numerical simulations without high-pass filtering show 
that these compact ripples present internal shock layers [20,24]. The main 
drawback of current (filtered) pseudospectral methods is the inability to show 
high-frequency phenomena. Particle methods based on the dispersive-velocity 
method have been proposed to cope with these features [22], but their preser- 
vation of the positivity of the solutions is another clear disadvantage, since 
after compacton collisions the solution may change sign. 

Both finite element and finite difference methods without high-frequency fil- 
tering have also been proposed. In finite element methods both a Petrov- 
Galerkin method using the product approximation developed by Sanz-Serna 
and Christie [20], and a standard method based on piecewise polynomials dis- 
continuous at the element interfaces [25] have been used. Second-order finite 
difference methods [26,27], high-order Pade methods [24], and the method of 
lines with adaptive mesh refinement [21] has also been applied with success. 
These methods also require artificial dissipation to simulate the generation of 
shocks after compacton interactions, which is usually incorporated by a lin- 
ear fourth-order derivative term. Such term introduces a plateau tail whose 
amplitude has been calculated for the K{2, 2) equation by Pikovsky and Rose- 
nau [15] by means of a variational perturbation theory for compactons. 

The main drawback in the numerical simulation of compacton propagation 
without high-frequency filtering is the appearance of spurious radiation, even 
in one-compacton solutions, as first shown by the authors in Ref. [28] by means 
of using the fourth-order Petrov-Galerkin finite element method developed by 
de Frutos, Lopez-Marcos and Sanz-Serna [20]. Both backward and forward 
propagating wavepackets of radiation are emitted from the compacton, having 
a very small amplitude, in fact, more than six orders of magnitude smaller than 
the compacton amplitude in current simulations. 

The main goal of this paper is a detailed analysis of the numerical origin 
and the main properties of the radiation emitted by compactons observed in 
Ref. [28]. First, in order to illustrate the universality of this phenomena, three 
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additional numerical methods are considered: the second order finite difference 
method developed by Ismail and Taha [26] and two Pade methods of sixth 
and eighth order developed by Rus and Villatoro [24]. Second, the numerical 
origin of the radiation is clarified by means of a grid refinement study. Third, 
in order to check if the origin of the radiation is due to the jump at the edge of 
the compacton suffered by its second-order derivative, the K{p,p) equations 
having compactons showing jumps in their first- to eight-order derivatives, i.e., 
with p G {3,2,5/3,3/2,7/5,4/3,9/7,5/4}, are also considered. And fourth, 
the graphical illustration of the self-similarity of the radiation is complemented 
with the numerical determination of their front velocity, wavepacket mean 
amplitude, and self-similar scaling exponents. 

The contents of this paper are as follows. Next Section presents the four nu- 
merical methods for the Rosenau-Hyman K{p,p) equation analyzed in this 
paper. Section 3 presents our results on the properties characterizing both the 
forward and the backward numerically-induced radiation wavepackets gen- 
erated during the propagation of one-compacton solutions. Finally, the last 
section is devoted to some conclusions. 



2 Numerical methods 



Let us consider the numerical solution of the RH Eq. (1) by means of the 
method of lines in time and several Pade approximations in space. Periodic 
boundary conditions in the interval x e [0,1/] are used as an approximation 
of the initial value problem in the whole real line. Let us take the fixed grid 
spacing Ax = L/M, the nodes Xm — rn Ax, for m = 0, 1, . . . , M, and a general 
Pade method written as 

A(E) ^ - coi3,(E) iUr.) + i3,(E) {UmY + C,(E) {Umf = 0, (3) 



where Um ~ u{xm), E is the shift operator, i.e., EUm = Um+i, the first and 
second derivatives are rationally approximated by means of Bi{E)/Ai{E) and 
Ci{E)/Ai(E), respectively, and i indicates the method among those studied in 
this paper. 

Method 1. The finite difference method developed by Ismail and Taha [26] 
is given by 



A(E)=T, 

_E-i + E^ 



Bi(E) 



2 Ax 



4 



^'^^^ = 2A^S ^ 



where X is the identity operator. In this case, method (3) is second-order 
accurate in space since 



Ai(E) dx 6 dx^ 

and 

Ci(E) Ax^ ff'u 
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Method 2. The finite element method developed by de Prutos et al. [20] is 
obtained by using 

, E-2 + 26E-^ + 66 + 26E^ + E2 
A(E) = 



B2(E) 



120 

-E-2 - 10 E"^ + lOE^ + E2 
24 Ax ' 



and C2(E) = Ci(E), where i32(E)/^2(E) and C2(E)/^2(E) are sixth- and 
fourth-order approximations to, respectively, the first- and third-order deriva- 
tives in Eq. (1), in fact 



.B2(E) du Ax^ d^u 

u — 



0(Ax«), 



A2{E) dx 5040 dx'^ 
and 

^u-^-^^+0(Ax'] 
A2{E) dx^ 240 dx''^^^ r 

Hence this method is fourth-order accurate in space. 

Method 3. A Pade method introduced in Ref. [24] which approximates the 
third- and first-order derivatives with, respectively, sixth- and fourth-order of 
accuracy, given by 

, E-2 + 56E-^ + 126 + 56E^ + E2 
•^"f^' = 240 ■ 
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B3(E) = B2(E), and C3(E) = Ci(E). In fact, Taylor series expansion yields 
^3(E) dx^ 240 9x5 + ^V^^^ 



and 

C,{E) d^u Ax' d% . . 
AsiE) dx^ 60480 dx^ ^ ^ 



Method 4. Another Pade method also introduced in Ref. [24] with an eighth- 
order accurate approximation to the first derivative in Eq. (1), obtained by 
means of 

, E-^ + 16E-^ + 36 + 16E^ + E2 
A(E) = ^ , 

^ -5 E-2 - 32 E-^ + 32 E^ + 5 E^ 
^'^^^ = 84A^ ' 



and C4(E) = Ci(E). This method is only of second-order for the third-order 
derivative, as shown by Taylor series expansion. Concretely, 



BJE) du Ax^ d^u 
A^iE) dx 44100 dx^ 
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u 



dx^ 28 dx^ 



0{Ax^) 



For sufficiently regular solutions of Eq. (1), Methods 2 and 3 arc fourth-order 
accurate, and Methods 1 and 4 only of second-order. Here on. Methods 1-4 
are referred to as Ismail, de Prutos, Pade-6, and Pade-8, respectively. Note 
that Methods 1-4 may be classified in function of the numerical order of 
approximation for the first and third derivatives in its local truncation error 
terms as (2,2), (6,4), (4,6), and (8,2), respectively. 

In this paper, the integration in time of Equation (3) is obtained by means of 
both the trapezoidal rule. 
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and the implicit midpoint rule, 
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where t"- = nAt and ~ u{xm, t""). Both methods are second-order accurate 
in time and yields implicit equations solved by using the Newton's method. 

The linear stability analysis by the von Neumann method for the methods 
developed in this section applied to the linearization of Eq. (1) shows its un- 
conditional (linear) stability [20,24,26]. Note that the usefulness of this linear 
stability analysis may be criticized when applied to a highly nonlinear prob- 
lem as Eq. (1), however, it is standard in a numerical analysis context. In fact, 
the solution of the four methods may blow-up for some Ax and At due to 
nonlinear instabilities whose analysis is outside the scope of this paper. 

Equation (1) has four invariants Ij = J (j)j{u)dx, where 0i = m, 02 = u^~^^, 
03 = M cos(a;), and (1)4 = u sin(a;). Methods 1-4 preserve exactly the first 
invariant of the K{p,p) equation, however, the other three invariants are not 
exactly preserved, but instead only well preserved [20,24,26]. 



3 Presentation of results 

Extensive numerical simulations of the K{p,p) equation with several p using 
either the trapezoidal or the implicit midpoint rule yield practically the same 
results for all of Methods 1-4, at least for At > Ax/10, hence, only results 
using the implicit midpoint rule are hereafter presented and discussed. For the 
sake of brevity, unless anything else is stated, the following figures and tables 
only show the results for the K{2, 2) equationPI 

The four plots in Fig. 1 show vertical zooms of the solution of the K{2, 2) 
equation at t = 300 for an initial condition given by one compacton with 
velocity c = 1 initially located at xq = 400, 500, 850, and 720 for, respectively, 
Ismail (top left plot), de Frutos (top right plot), Pade-6 (bottom left plot), 
and Pade-8 (bottom right plot) methods. The four plots in Fig. 1 clearly 
show that two wavepackets of radiation are generated from the compacton, 
here referred to as forward and backward radiation corresponding to that 

^ Supplementary material with figures and tables presenting results for the K{p,p) 
equation with p € {3,2,5/3,3/2,7/5,4/3,9/7,5/4} may be found in the web page 
http : / / www . Icc . uma . es/ ~rusmaii/iiivest/compact/Compactons . htm. 
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Fig. 1. Zoom in of snapshots at t = 300 of the radiation generated at both sides of 
a compacton of the K{2, 2) equation propagating with c = 1 numerically calculated 
using Co = 1, Ax = 0.05, and At = 0.1 by means of Ismail (top left plot), de Frutos 
(top right plot), Pade-6 (bottom left plot), and Pade-8 (bottom right plot) methods, 
initially located at, respectively, xo = 400, 500, 850, and 720. 



propagating to the right and to the left, respectively, of the compacton. Note 
that the initial position of the compactons is not the same in all the plots 
in order to avoid that the backward (forward) radiation cross the left (right) 
boundary reappearing through the other one due to the periodic boundary 
conditions used in the simulations. Note also the use of cq = c in order to 
stop the compacton and highhght the relative velocity of both wavepackets of 
radiation generated during its propagation. 

The plots in Fig. 1 show that the amplitude of both wavepackets is very small 
compared with that of the compacton, being that of the backward radiation 
two orders of magnitude larger than that of the forward one for Ismail (Fig. 1, 
top left plot) and Pade-8 (bottom right plot), but only several times largest for 
the other two methods. The backward radiation has a steeper front than that 
of the forward one for all the methods and a front velocity smaller (in absolute 
value) than the forward one for de Frutos (top right plot), Pade-6 (bottom 
left plot), and Pade-8 (bottom right plot) methods, being approximately equal 
for Ismail (top left plot) one. Note that, for long time integrations under pe- 
riodic boundary conditions, both radiation wavepackets collide resulting in a 
background dominated by the backward radiation, whereon the compacton 
propagates, due to its robustness, without appreciable change on its parame- 
ters. 
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The more interesting and noticeable property of both backward and forward 
compacton radiation is their self-similarity. Figures 2, 3, 4, and 5 show the 
absolute value of both the forward (right plots) and the backward (left ones) 
radiation for, respectively, Ismail, de Prutos, Pade-6, and Pade-8 methods at 
time t — 150 (top plots) and t — 300 (bottom ones). In the plots of Figs. 2-5, 
the horizontal axis is selected in order to best illustrate the self- similarity of 
the wavepacket envelope of both numerically-induced radiations by graphical 
comparison of the top and bottom plots. As shown in Figs. 2-5 the envelope 
shape of both the forward and the backward radiations is highly dependent 
on the method and, not illustrated in the plots, on its parameters Ax and At, 
and the compacton velocity c. 

A possible origin of the self-similar radiation may be the jump experienced 
by the second-order derivative of the K{2,2) compacton at their edges. How- 
ever, extensive numerical simulations using Methods 1-4 show that the self- 
similarity of the radiation is also present in the propagation of K{p,p) com- 
pactons with k continuous derivatives at its both edges, i.e., for p = {2 + k)/k. 
The only cases in which the self-similarity is not clearly visible are for p ~ 1 , 
for which the amplitude of the radiation is comparable with the tolerance used 
in the iterations of the Newton method. In such cases, the radiation near the 
compacton is degraded by noise, apparently introduced by round-off errors, 
destroying the self-similarity and, for long-time integrations, blowing up the 
solution. 

The numerical origin of the spurious radiation observed in the simulations is 
illustrated in Tables 1 and 2, which show the amplitude of both the backward 
(ub) and forward (uj) radiation for the four methods studied in this paper as 
a function of Ax and At, respectively. This amplitude has been determined 
by finding the first local maximum of the wavepacket starting from the front 
of the wavepacket using a five point rule, i.e., three nodes where the function 
increases followed by two nodes where it decreases, being the amplitude value 
that of the central node. Tables 1 and 2 clearly show that the amplitude of 
both radiations decreases with decreasing Ax but remains practically constant 
with At. The numerical origin of the radiations appear to be the numerical 
approximation of the spatial derivatives. The last column of Table 1 shows 
the exponent q such that the amplitude of the radiations are 0{Ax'^), calcu- 
lated by means of linear regression. This exponent may clarified whether the 
approximation of either the first or the third derivatives in Eq. (1) is the only 
responsible of any of these radiations. However, the results shown in Table 1 
are not conclusive in this respect and the radiations appear to be the result 
of the trade-off between the local truncation error of both derivatives. Simi- 
lar results have been obtained for the other K{p,p) equations studied in this 
paper. 

The position of the left (right) front of the backward (forward) radiation 
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Ismail: 1-c. back. rad. Ismail: 1-c. tor. rad. 
' ' 1 10"'f ■ ' 




Fig. 2. Backward (left plots) and forward (right ones) radiation generated by a 
compacton of the ^(2,2) equation numerically propagating with Ismail method 
with Ax = 0.05, At = 0.1, and cq = c = 1 at two instants of time, t = 150 (top 
plots) and t = 300 (bottom ones), highlighting their self-similarity. 

Table 1 

Amplitude at t = 150 of both the backward (tif,) and forward (uf) radiation for a 
compacton of the K{2, 2) equation with velocity c = 1 as a function of Ax using 
At = 0.05, Co = c, and x G [0,2500]. The asterisks indicate solutions which blow 
up. Linear regression is used to obtain q such that Uf and Uf, are O(Ax^). 



Meth. 
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0.2 


0.1 


0.05 


0.025 


0.0125 
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Uf 
Ub 


6.54 X 10"^ 
2.78 X 10-^ 


7.14 X 10-6 
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5.68 X 10-6 
1.45 X 10-^ 
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2.11 X 10-^ 
7.43 X 10-^ 


8.12 X 10-9 
2.35 X 10-^ 


2.4 
2.5 
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Ub 


1.76 X 10"^ 
4.24 X 10-3 


1.64 X 10^6 
3.17 X lO-'' 


2.79 X lO-'^ 
8.90 X 10-5 


4.84 X 10"^ 
2.34 X 10-5 


7.23 X 10-9 
5.97 X 10-6 
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de Frutos: 1-c. back. rad. , deFrutos: 1-c. for. rad. 




Fig. 3. Backward (left plots) and forward (right ones) radiation generated by a 
compacton of the K{2, 2) equation numerically propagating with de Prutos method 
with Ax = 0.05, At = 0.1, and cq = c = 1 at two instants of time, t = 150 (top 
plots) and t = 300 (bottom ones), highlighting their self-similarity. 

Table 2 

Amplitude at i = 150 of both the backward («(,) and forward {uf) radiation for a 
compacton of the K{2, 2) equation with velocity c = 1 as a function of At using 
Co = c. Ax = 0.05, and x G [0, 2500]. Linear regression is used to obtain q such that 
Uf and Ub are 0{At'^). 



Meth. 


At 


0.1 


0.05 


0.025 


0.0125 


0.00625 
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Uf 


1.04 X 10-6 
2.17 X 10-^ 


1.26 X 10-6 
2.94 X 10-^ 


1.37 X 10-6 
1.53 X 10-^ 


1.40 X 10-6 
1.44 X 10-3 


1.41 X 10-6 
1.21 X 10-3 
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-.73 
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1.52 X 10-^ 
2.72 X 10-^ 


1.61 X 10-^ 
2.60 X 10-^ 


1.96 X 10-^ 
4.18 X 10-^ 


2.11 X 10-^ 
7.88 X 10-'^ 


2.16 X 10-'^ 
7.70 X 10-'^ 


-.14 
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1.53 X 10-^ 
4.11 X 10-^ 


8.95 X 10-8 
3.81 X 10-'^ 


1.11 X 10-^ 
1.11 X 10-6 


1.24 X 10-^ 
8.71 X 10-^ 


1.29 X 10-'^ 
7.00 X 10-'^ 


.002 
-.27 


4 


Ub 


2.65 X lO"'^ 
8.64 X 10"^ 


2.79 X lO"'^ 
8.90 X 10-^ 


2.75 X 10""^ 
8.18 X 10-^ 


2.81 X lO-'^ 
8.17 X 10-^ 


2.93 X lO-'^ 
8.08 X 10-^ 


-.03 
.030 
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Pade-6: 1 -c. back. rad. , Pade-6: 1 -c. for. rad. 
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Fig. 4. Backward (left plots) and forward (right ones) radiation generated by a 
compacton of the K{2, 2) equation numerically propagating with Pade-6 method 
with Ax = 0.05, At = 0.1, and cq = c = 1 at two instants of time, t = 150 (top 
plots) and t = 300 (bottom ones), highlighting their self-similarity. 

wavepackets for the Ismail (top left plot), de Prutos (top right plot), Pade-6 
(bottom left plot) and Padc-8 (bottom right plot) methods for the K{2, 2) 
equation is shown in Fig. 6. This position has been determined, using linear 
interpolation, as the "first" point from the outside of the wavepacket, i.e., 
from left to right (right to left) for backward (forward) radiation, where the 
amplitude of the solution is equal to an amplitude threshold, the half of the 
maximum amplitude of the radiation at t = 300. The four plots in Fig. 6 
clearly show a linear evolution of the position of the front for both forward 
(continuous line) and backward (dashed line) radiations. The velocity of the 
front, i.e., the slope of these curves, is nearly constant during propagation 
being negative (positive) for backward (forward) radiation. The constancy of 
the front velocities has also been observed in the simulations of the K{p,p) 
equation. 

The front velocity of both the forward (c/) and the backward (c^) wavepackets 
relative to the velocity of the compacton may be calculated by linear regression 
from the evolution in time of their positions. Our extensive numerical exper- 
iments show that both the forward (cy) and backward (cf,) front velocities 
are linear functions of the parameter Cq, being practically independent of the 
parameters Ax and At, and, also nearly independent of c and p, for a K{p,p) 
compacton, with a small percentage increase as p decreases approaching unity. 
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Pade-8: 1 -c. back. rad. 




Pade-8: 1 -c. for. rad. 
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Fig. 5. Backward (left plots) and forward (right ones) radiation generated by a 
compacton of the K{2, 2) equation numerically propagating with Pade-8 method 
with Ax = 0.05, At = 0.1, and cq = c = 1 at two instants of time, t = 150 (top 
plots) and t = 300 (bottom ones), highlighting their self-similarity. 
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Fig. 6. Front velocity evolution of both forward (continuous line) and backward 
(dashed line) radiations for a K{2, 2) compacton numerically propagating as a func- 
tion of time for Ismail (top left plot), de Frutos (top right), Pade-6 (bottom left) 
and Pade-8 (bottom right) methods with Ax = 0.05, At = 0.1 and cq = c = 1. 
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Table 3 

Velocity of the front of both the backward and forward radiation for a K{2, 2) 
compacton with velocity c = 1 as function of cq in numerical simulations with two 
Ax, two At, X e [0, 2500] and t G [0, 100]. 
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1/2 


1 


2 


1/2 


1 


2 


1/2 


1 


2 


1 




* 


1.02 


2.03 


* 


1.02 


2.03 


* 


1.01 


2.02 




Cb 


* 


-1.02 


-2.03 


* 


-1.03 


-2.04 


* 


-1.02 


-2.03 


2 




2.55 


5.06 


10.1 


2.56 


5.07 


10.1 


2.53 


5.05 


* 




Cb 


-0.500 


-1.01 


-2.01 


-0.504 


-1.01 


-2.02 


-0.505 


-1.01 


* 


3 


c/ 


5.09 


10.1 


20.1 


5.09 


10.1 


20.2 


5.06 


10.1 


* 




Cb 


-0.672 


-1.25 


-2.05 


-0.634 


-1.03 


-2.02 


-0.514 


-1.01 


* 


4 


c/ 


3.12 


6.18 


12.3 


3.12 


6.19 


12.4 


3.10 


6.16 


* 




Cb 


-0.510 


-1.00 


-2.01 


-0.506 


-1.01 


-2.02 


-0.505 


-1.01 


* 



Table 3 shows both the forward (cf) and backward (c{,) front velocities for 
Co = c/2, c, and 2c, for two values of Ax and two values of At. The front 
velocities depend linearly on Cq instead on c, in fact, Cf ~ Cq, Scq, 10 Cq, 
and 6 Cq for Methods 1-4, respectively, and c^ ~ — Cq for the four methods. 
This approximations are better as Ax decreases. Let us highlight that both 
front velocities are relative to that of the compacton, therefore, in a rest frame 
of reference, where the compacton propagates with its own velocity instead of 
being stopped by the condition Cq = c, the backward radiation is generated 
in the left edge of the compacton at t = and stretches as it propagates, like 
a wake left in the track of the compacton during its propagation signaling its 
initial position in the numerical simulation. 

Present results suggest that the spatial numerical approximation of the lin- 
ear term introduced in Eq. (1) in order to stop de compacton may be the 
responsible of the self-similarity of the envelope of the radiation studied in 
this paper. Introducing into the modified equation [29] for Method i applied 
to Eq. (1) the solution u = Uc + Ur, where Uc is the compacton (2) and Ur 
is the numerically induced radiation, with \ur\ -C \uc\ in the support of the 
compacton and Uc = outside it, yields 
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Fig. 7. Group velocity of Method 1 (solid line), 2 (dashed line), 3 (dash-dotted line), 
and 4 (dotted line) for the linear Eq. (6) as a function of the normalized wavenumber 
a, such that k = air /Ax. 

where D = d/dx and h.o.t. stands for higher-order terms. The hnear disper- 
sion of this equation is obtained by substitution of = exp(i {kx — Wi{k) t)) 
into Eq. (6). The envelope of a wavepacket of radiation propagates with the 
group velocity, Ci{k) = dwi{k)/dk, given by 



„ d sm{/\xk) , . , , 

Ci{k) = -Co — — ^ — = -Co cos(AxA;), 

^ d 50 sin(Aa; /c) + 5 sin(2 Ax A;) 

' ~ d^(33 + 26 cos(Aa;A;) + cos(2AxA;)) Ax' 

^ ^ d 100 sin(Aa; A;) + 10 sin(2 Ax A;) 

'~~^dk (63 + 56 cos( Ax k) + cos(2 Ax A;)) Ax ' 

d 160 sin(Ax k) + 25 sin(2 Ax k) 
' ~~^^dk (108 + 96 cos(Ax k)+6 cos(2 Ax k)) Ax' 

for Methods 1-4, respectively, which are plotted in Figure 7 as a function 
of the normalized wavenumber a, given by A; = a k^^ax where the highest 
wavenumber in the spatial grid is A;max = vr/Ax. The discrete Fourier transform 
of the forward radiation shows that its spectrum is concentrated around the 
highest wavenumber, therefore, its front velocity is given by ^(/Cinax) = Cq, 
5 Co, 10 Cq, and 6.11 Co, for Methods 1-4, respectively. The spectrum for the 
backward radiation presents several peaks of low frequency accompanied with 
a smaller peak at the highest wavenumber, therefore, its front velocity is given 
approximately by Ci(/cmax/10) ~ — cq, for Methods 1-4. These results are in 
good agreement with Table 3 and further results for the K{p,p) equation 
omitted here for brevity. 
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Equation (6) is not the only responsible of the numerically induced self-similar 
radiations by the K{p,p) compactons since it is easy to show that it has no 
self-similar solutions. Moreover, the numerical solution by all the methods 
studied in this paper blow ups when cq either has a negative value (cq c < 0) 
or has a value very different from c (either |co| <S |c|, or |co| ^ |c|). This result, 
found in a large number of simulations and illustrated in the last cohimn of 
Table 3, was unnoticed by the authors of Ref. [20], whose first introduced 
the linear term in Eq. (1), and in further references [24,26,27]. Furthermore, 
a linear stability analysis of the semidiscrete Eq. (6) shows its unconditional 
stability, independently of the value of Cq, so the roots of this instability must 
be in the nonlinear terms not considered in it. 

The self-similarity of both the forward and the backward radiation requires 
that their analytical expressions be self-similar functions which may be ana- 
lytically written as, respectively. 



where Xb = Xo — tt / {2 (3) and x f = Xo + / {2 f3) are, respectively, the left and 
right extremes of the compacton solution, Qb and Qf are the scaling exponents 
for, respectively, the forward and the backward radiation, and Uf and Ub are 
the shapes of, respectively, the forward and backward wavepackets. In order 
to verify that the scaling exponents in Eqs. (7) and (8) are really constant, 
the temporal evolution of the amplitude of both radiations must be studied. 
Figure 8 shows that this amplitude changes a little in the first steps of time but 
yields a very smooth decreasing curve as time progresses, being approximately 
linear in the logarithmic scale of the plots. Therefore, the temporal evolution 
of the amplitude is asymptotically exponential in time. Similar results have 
been also obtained for other K{p,p) compactons and/or mesh parameters. 

Tables 4 and 5 show the scahng exponents in Eqs. (7) and (8) as function 
of Ax and Cq = c, respectively, since our extensive set of simulations show 
that the time step has no significant influence on the results. The scaling 
exponents in these tables have been determined by using linear regression 
of the temporal evolution of the "mean" amplitude of the envelope of the 
wavepackets, i.e., the mean of the absolute value of the amplitude of the 
solution in the intervals [xb,Xb + ct] and [xf + ct,Xf + {c + Cf)t] for the 
backward and forward radiations, respectively. To avoid the effects of the 
initial transient, where the self-similarity of the wavepackets is not properly 
defined due to the aliasing errors introduced by the sampling of the solution, 
the first 25% of the solution is not considered in the linear regression. 




— oo < X < Xb + ct, 



Xf + ct < X < oo. 



(7) 



(8) 
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Ismail de Frutos 




100 200 300 100 200 

t t 



Fig. 8. Amplitude evolution of both forward (continuous line) and backward (dashed 
line) radiations as a function of time for Ismail (top left plot), de Frutos (top right), 
Padc-6 (bottom left) and Pade-8 (bottom right) methods applied to the K{2,2) 
equation. 



Table 4 

Scaling exponent of both forward and backward radiations for the Ismail, de Frutos, 
Pade-6, and Pade-8 methods for a K{2, 2) compacton in numerical simulations with 
At = 0.05 and cq = c = 1 as a function of Ax, calculated using a linear regression 
of the evolution in time of the mean amplitude in the interval t € [75, 300]. 





Ax 


0.2 


0.1 


0.05 


0.025 


0.0125 


Ismail 




0.487 


0.453 


0.503 


0.482 


0.464 




Qb 


0.487 


0.924 


0.977 


1.05 


1.18 


de Frutos 


Bf 


0.501 


0.497 


0.498 


0.499 


0.500 




Qh 


0.549 


0.547 


0.514 


0.513 


0.505 


Pade-6 


Bf 


0.492 


0.494 


0.496 


0.498 


0.608 




Qh 


0.532 


0.537 


0.500 


0.506 


0.612 


Pade-8 


Bf 


0.481 


0.513 


0.527 


0.556 


0.520 




Qh 


0.077 


0.774 


0.931 


0.891 


0.875 
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Table 5 

Scaling exponent of both forward and backward radiations for the Ismail, de Frutos, 
Pade-6, and Pade-8 methods for a K{2, 2) compacton in numerical simulations with 
At = 0.05 and Ax = 0.05 as a function of cq = c, calculated using a linear regression 
of the evolution in time of the mean amplitude in the interval t G [75, 300]. 





Co 


0.1 


0.2 


0.5 


1 


1.5 


2 


5 


Ismail 




0.157 


0.512 


0.519 


0.503 


0.492 


0.487 


0.477 




Qb 


2.408 


1.18 


1.07 


0.977 


0.827 


0.794 


0.778 


de Frutos 


Sf 


0.492 


0.494 


0.496 


0.510 


0.495 


0.495 


0.498 




Qb 


0.542 


0.552 


0.505 


0.511 


0.499 


0.495 


0.479 


Pade-6 


Sf 


0.487 


0.489 


0.491 


0.496 


0.489 


0.490 


0.489 




Qb 


0.561 


0.507 


0.437 


0.500 


0.496 


0.507 


0.504 


Pade-8 


Qf 


0.534 


0.496 


0.562 


0.527 


0.530 


0.507 


0.482 




Qb 


0.433 


0.962 


0.995 


0.931 


0.899 


0.760 


0.483 



Tables 4 and 5 show that the scaling factor Qf is approximately equal to 0.5 
for all the four methods studied in this paper, nearly independent of both 
Ax and cq = c, respectively, with the largest dispersion associated to Pade-8 
and Ismail methods. Tables 4 and 5 also show that the scahng factor is 
approximately equal to 0.5 for the de Frutos and Pade-6 methods, to 0.9 for 
the Pade-8 method, and to 1.0 for the Ismail method, with a small decrement 
as c grows, for both the Pade-8 and Ismail methods. For the last two methods 
the dispersion is large, although diminish if the first column of Tables 4 and 5 
is not taken into account since it corresponds to a large Ax and presents very 
noticeable aliasing errors. 



4 Conclusions 

The propagation of compactons of the Rosenau-Hyman K{p,p) equation has 
been studied by means of four numerical methods showing the appearance of 
numerically-induced radiation. Both backward and forward wavcpackcts are 
generated from the compacton with a clear self- similar shape, illustrated by 
means of properly scaling of the figures. The parameters characterizing these 
wavepackets have been numerically determined. An analytical model, based 
on the linearization of K{p,p) equation has been used in order to approxi- 
mate the front velocities of the forward wavcpackcts relative to that of the 
compacton. The front velocity of the backward wavepacket is approximately 
equal to that of the compacton but with opposite sign. Both forward and 
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backward front velocities are practically independent of the parameters Ax 
and At of the numerical method. The evolution in time of the mean ampli- 
tude of the wavepackets shows its exponential decreasing in time, suggesting 
an ansatz for both the backward and forward wavepackets corresponding to 
self-similar functions characterized by scaling exponents, which have been nu- 
merically calculated for both radiations by means of the linear regression of 
the logarithm of the mean amplitude, showing that its value is approximately 
constant as a function of both the mesh grid size and the compacton velocity. 

The scaling exponents approximately equal to 1/2 for the amplitude evolution 
of the envelope of both the forward and the backward radiations suggest that 
they may be analytically approximated for weak nonlinearity (valid for the 
radiation but not for the compactons) by a nonlinear (cubic) Schrodinger 
equation. This analysis is in progress. Asymptotic analysis using the method 
of modified equations applied to the four numerical methods studied in this 
paper may also be useful in the characterization of the radiation wavepackets 
found here. In fact, the analytical explanation of both the generation of the 
self-similar radiations and the blow-up of the solution for some Cq are very 
interesting open problems. 
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